Предобработка данных

In [1]:
# Подключение к Google drive

from google.colab import drive
drive.mount('/content/drive')
Go to this URL in a browser: https://accounts.google.com/o/oauth2/auth?client_id=947318989803-6bn6qk8qdgf4n4g3pfee6491hc0brc4i.apps.googleusercontent.com&redirect_uri=urn%3aietf%3awg%3aoauth%3a2.0%3aoob&scope=email%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdocs.test%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive%20https%3a%2f%2fwww.googleapis.com%2fauth%2fdrive.photos.readonly%20https%3a%2f%2fwww.googleapis.com%2fauth%2fpeopleapi.readonly&response_type=code

Enter your authorization code:
··········
Mounted at /content/drive
In [2]:
import numpy as np
import pandas as pd

Специфические преобразования

In [3]:
# Загрузим набор данных

df = pd.read_csv('/content/drive/My Drive/Colab Notebooks/freMPL-R.csv', low_memory=False)
df = df.loc[df.Dataset.isin([5, 6, 7, 8, 9])]
df.drop('Dataset', axis=1, inplace=True)
df.dropna(axis=1, how='all', inplace=True)
df.drop_duplicates(inplace=True)
df.reset_index(drop=True, inplace=True)
In [3]:
df = pd.read_parquet('data/df.parquet.gzip', engine = 'fastparquet')

В предыдущем уроке мы заметили отрицательную величину убытка для некоторых наблюдений. Заметим, что для всех таких полисов переменная "ClaimInd" принимает только значение 0. Поэтому заменим все соответствующие значения "ClaimAmount" нулями.

In [4]:
NegClaimAmount = df.loc[df.ClaimAmount < 0, ['ClaimAmount','ClaimInd']]
print('Unique values of ClaimInd:', NegClaimAmount.ClaimInd.unique())
NegClaimAmount.head()
Unique values of ClaimInd: [0]
Out[4]:
ClaimAmount ClaimInd
index
145862 -74.206042 0
145955 -1222.585196 0
145957 -316.288822 0
146143 -666.758610 0
146155 -1201.600604 0
In [5]:
df.loc[df.ClaimAmount < 0, 'ClaimAmount'] = 0

Для моделирования частоты убытков сгенерируем показатель как сумму индикатора того, что убыток произошел ("ClaimInd") и количества заявленных убытков по различным видам ущерба за 4 предшествующих года ("ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen").

В случаях, если соответствующая величина убытка равняется нулю, сгенерированную частоту также обнулим.

In [6]:
df['ClaimsCount'] = df.ClaimInd + df.ClaimNbResp + df.ClaimNbNonResp + df.ClaimNbParking + df.ClaimNbFireTheft + df.ClaimNbWindscreen
df.loc[df.ClaimAmount == 0, 'ClaimsCount'] = 0
df.drop(["ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen"], axis=1, inplace=True)
In [7]:
pd.DataFrame(df.ClaimsCount.value_counts()).rename({'ClaimsCount': 'Policies'}, axis=1)
Out[7]:
Policies
0.0 104286
2.0 3529
1.0 3339
3.0 2310
4.0 1101
5.0 428
6.0 127
7.0 26
8.0 6
9.0 2
11.0 1
In [8]:
import plotly.express as px
fig = px.scatter(df, x='ClaimsCount', y='ClaimAmount', title='Зависимость между частотой и величиной убытков')
fig.show()

Для моделирования среднего убытка можем рассчитать его как отношение величины убытков к их частоте.

In [9]:
df.loc[df.ClaimsCount > 0, 'AvgClaim'] = df['ClaimAmount']/df['ClaimsCount']

Общие преобразования

Класс для общих случаев

In [10]:
class InsDataFrame:


    ''' Load data method '''

    def load_pd(self, pd_dataframe):
        self._df = pd_dataframe


    ''' Columns match method '''

    def columns_match(self, match_from_to):
        self._df.rename(columns=match_from_to, inplace=True)


    ''' Person data methods '''

    # Gender
    _gender_dict = {'Male':0, 'Female':1}

    def transform_gender(self):
        self._df['Gender'] = self._df['Gender'].map(self._gender_dict)

        

    # Age
    @staticmethod
    def _age(age, age_max):
        if pd.isnull(age):
            age = None
        elif age < 18:
            age = None
        elif age > age_max:
            age = age_max
        return age
      
    def transform_age(self, age_max=70):
        self._df['driver_minage'] = self._df['driver_minage'].apply(self._age, args=(age_max,))

    # Age M/F
    @staticmethod
    def _age_gender(age_gender):
        _age = age_gender[0]
        _gender = age_gender[1]
        if _gender == 0: #Male
            _driver_minage_m = _age
            _driver_minage_f = 18
        elif _gender == 1: #Female
            _driver_minage_m = 18
            _driver_minage_f = _age
        else:
            _driver_minage_m = 18
            _driver_minage_f = 18
        return [_driver_minage_m, _driver_minage_f]
    
    def transform_age_gender(self):
        self._df['driver_minage_m'],self._df['driver_minage_f'] = zip(*self._df[['driver_minage','Gender']].apply(self._age_gender, axis=1).to_frame()[0])

    # Experience
    @staticmethod
    def _exp(exp, exp_max):
        if pd.isnull(exp):
            exp = None
        elif exp < 0:
            exp = None
        elif exp > exp_max:
            exp = exp_max
        return exp

    def transform_exp(self, exp_max=52):
        self._df['driver_minexp'] = self._df['driver_minexp'].apply(self._exp, args=(exp_max,))


    ''' Other data methods '''

    def polynomizer(self, column, n=2):
        if column in list(self._df.columns):
            for i in range(2,n+1):
                self._df[column+'_'+str(i)] = self._df[column]**i

    def get_dummies(self, columns):
        self._df = pd.get_dummies(self._df, columns=columns)


    ''' General methods '''

    def info(self):
        return self._df.info()

    def head(self, columns, n=5):
        return self._df.head(n)

    def len(self):
        return len(self._df)

    def get_pd(self, columns):
        return self._df[columns]

Создаем класс-наследник, в котором переопределяем некоторые методы, специфические для конкретной ситуации, и создаем новые

  • В данных стаж "LicAge" измеряется в неделях.
  • Фактор "SocioCateg" содержит информацию о социальной категории в виде кодов классификации CSP. Агрегируем имеющиеся коды до 1 знака, а затем закодируем их с помощью one-hot encoding.

Wiki

Более подробный классификатор

In [11]:
class InsDataFrame_Fr(InsDataFrame):

    # Experience (weeks to years)
    @staticmethod
    def _exp(exp, exp_max):
        if pd.isnull(exp):
            exp = None
        elif exp < 0:
            exp = None
        else:
            exp * 7 // 365
        if exp > exp_max:
            exp = exp_max
        return exp

    # Marital status
    _MariStat_dict = {'Other':0, 'Alone':1}

    def transform_MariStat(self):
        self._df['MariStat'] = self._df['MariStat'].map(self._MariStat_dict)
    
    # Social category
    def transform_SocioCateg(self):
        self._df['SocioCateg'] = self._df['SocioCateg'].str.slice(0,4)

Подгружаем данные

In [12]:
data = InsDataFrame_Fr()
In [13]:
data.load_pd(df)
In [14]:
data.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 115155 entries, 145780 to 310979
Data columns (total 20 columns):
 #   Column        Non-Null Count   Dtype  
---  ------        --------------   -----  
 0   Exposure      115155 non-null  float64
 1   LicAge        115155 non-null  int64  
 2   RecordBeg     115155 non-null  object 
 3   RecordEnd     59455 non-null   object 
 4   Gender        115155 non-null  object 
 5   MariStat      115155 non-null  object 
 6   SocioCateg    115155 non-null  object 
 7   VehUsage      115155 non-null  object 
 8   DrivAge       115155 non-null  int64  
 9   HasKmLimit    115155 non-null  int64  
 10  BonusMalus    115155 non-null  int64  
 11  ClaimAmount   115155 non-null  float64
 12  ClaimInd      115155 non-null  int64  
 13  OutUseNb      115155 non-null  float64
 14  RiskArea      115155 non-null  float64
 15  PolicyCount   115155 non-null  int64  
 16  ClaimCount    115155 non-null  bool   
 17  NoClaimCount  115155 non-null  int64  
 18  ClaimsCount   115155 non-null  float64
 19  AvgClaim      10869 non-null   float64
dtypes: bool(1), float64(6), int64(7), object(6)
memory usage: 22.7+ MB

Преобразовываем параметры

In [15]:
# Переименовываем
data.columns_match({'DrivAge':'driver_minage','LicAge':'driver_minexp'})
In [16]:
# Преобразовываем
data.transform_age()
data.transform_exp()
data.transform_gender()
data.transform_MariStat()
data.transform_SocioCateg()
In [17]:
# Пересечение пола и возраста, их квадраты
data.transform_age_gender()
data.polynomizer('driver_minage_m')
data.polynomizer('driver_minage_f')

Для переменных, содержащих более 2 значений:

  • либо упорядочиваем значения,
  • либо используем фиктивные переменные (one-hot encoding).

NB: В H2O не рекомендуется использовать one-hot encoding. Тем не менее, используем здесь фиктивные переменные, чтобы в дальнейшем сохранить возможность сравнения результатов построенных моделей.

In [18]:
# Onehot encoding
data.get_dummies(['VehUsage','SocioCateg'])
In [19]:
data.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 115155 entries, 145780 to 310979
Data columns (total 33 columns):
 #   Column                           Non-Null Count   Dtype  
---  ------                           --------------   -----  
 0   Exposure                         115155 non-null  float64
 1   driver_minexp                    115155 non-null  int64  
 2   RecordBeg                        115155 non-null  object 
 3   RecordEnd                        59455 non-null   object 
 4   Gender                           115155 non-null  int64  
 5   MariStat                         115155 non-null  int64  
 6   driver_minage                    115155 non-null  int64  
 7   HasKmLimit                       115155 non-null  int64  
 8   BonusMalus                       115155 non-null  int64  
 9   ClaimAmount                      115155 non-null  float64
 10  ClaimInd                         115155 non-null  int64  
 11  OutUseNb                         115155 non-null  float64
 12  RiskArea                         115155 non-null  float64
 13  PolicyCount                      115155 non-null  int64  
 14  ClaimCount                       115155 non-null  bool   
 15  NoClaimCount                     115155 non-null  int64  
 16  ClaimsCount                      115155 non-null  float64
 17  AvgClaim                         10869 non-null   float64
 18  driver_minage_m                  115155 non-null  int64  
 19  driver_minage_f                  115155 non-null  int64  
 20  driver_minage_m_2                115155 non-null  int64  
 21  driver_minage_f_2                115155 non-null  int64  
 22  VehUsage_Private                 115155 non-null  uint8  
 23  VehUsage_Private+trip to office  115155 non-null  uint8  
 24  VehUsage_Professional            115155 non-null  uint8  
 25  VehUsage_Professional run        115155 non-null  uint8  
 26  SocioCateg_CSP1                  115155 non-null  uint8  
 27  SocioCateg_CSP2                  115155 non-null  uint8  
 28  SocioCateg_CSP3                  115155 non-null  uint8  
 29  SocioCateg_CSP4                  115155 non-null  uint8  
 30  SocioCateg_CSP5                  115155 non-null  uint8  
 31  SocioCateg_CSP6                  115155 non-null  uint8  
 32  SocioCateg_CSP7                  115155 non-null  uint8  
dtypes: bool(1), float64(6), int64(13), object(2), uint8(11)
memory usage: 25.6+ MB
In [20]:
col_features = [
                'driver_minexp',
                'Gender',
                'MariStat',
                'HasKmLimit',
                'BonusMalus',
                'OutUseNb',
                'RiskArea',
                'driver_minage_m',
                'driver_minage_f',
                'driver_minage_m_2',
                'driver_minage_f_2',
                'VehUsage_Private',
                'VehUsage_Private+trip to office',
                'VehUsage_Professional',
                'VehUsage_Professional run',
                'SocioCateg_CSP1',
                'SocioCateg_CSP2',
                'SocioCateg_CSP3',
                'SocioCateg_CSP4',
                'SocioCateg_CSP5',
                'SocioCateg_CSP6',
                'SocioCateg_CSP7'  
]
In [21]:
col_target = ['ClaimAmount', 'ClaimsCount', 'AvgClaim']
In [22]:
df_freq = data.get_pd(col_features+col_target)
In [23]:
df_ac = df_freq[df_freq['ClaimsCount'] > 0].reset_index().copy()

Разделение набора данных на обучающую, валидационную и тестовую выборки

In [24]:
from sklearn.model_selection import train_test_split
In [25]:
# Разбиение датасета для частоты на train/val/test

x_train_c, x_test_c, y_train_c, y_test_c = train_test_split(df_freq[col_features], df_freq.ClaimsCount, test_size=0.3, random_state=1)
x_valid_c, x_test_c, y_valid_c, y_test_c = train_test_split(x_test_c, y_test_c, test_size=0.5, random_state=1)
In [26]:
# Разбиение датасета для среднего убытка на train/val/test 

x_train_ac, x_test_ac, y_train_ac, y_test_ac = train_test_split(df_ac[col_features], df_ac.AvgClaim, test_size=0.3, random_state=1)
x_valid_ac, x_test_ac, y_valid_ac, y_test_ac = train_test_split(x_test_ac, y_test_ac, test_size=0.5, random_state=1)

Обобщенные линейные модели (Generalized Linear Models, GLM)

Теория

Пусть $y$ – целевая переменная, $X$ – матрица объясняющих переменных, $\beta$ – вектор параметров модели.

Матрица $X$ составлена из всех векторов наблюдений $x_i$, каждый из которых представляет собой объясняющую переменную.

Основные компоненты обобщенной линейной модели:

  • Систематическая компонента $\eta$:
    • $\eta = X\beta \hspace{10pt}(=\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_nx_n)$.
  • Случайная компонента $y$:
    • Элементы вектора $y$ – независимые одинаково распределенные случайные величины, имеющие функцию плотности распределения $f(y;\theta,\phi)$ из экспоненциального семейства.
    • Распределения из экспоненциального семейства имеют параметры $\theta$ (характеристика среднего) и $\phi$ (характеристика дисперсии). В общем виде данные распределения могут быть определены: $$f_i(y_i;\theta_i,\phi)=\exp\left\lbrace \frac{y_i\theta_i-b(\theta_i)}{a_i(\phi)} + c(y_i, \phi) \right\rbrace,$$ где $a_i(\phi)$, $b(\theta_i)$ и $c(y_i, \phi)$ некоторые функции.
    • Для распределений из данного семейства дисперсия является функцией от среднего.
    • Экспоненциальное семейство включает распределения нормальное, экспоненциальное, Пуассона, гамма, хи-квадрат, бета и другие.
  • Функция связи $g$:
    • $\mathbb{E}\left[y\right]=\mu=g^{-1}\left(\eta\right)$, $\mu$ – математическое ожидание $y$;
    • $g$ – монотонная дифференцируемая функция.

GLM с распределением Пуассона

Регрессия Пуассона обычно используется в случаях, когда зависимая переменная представляет собой счетные значения и ошибки предполагаются распределенными в соответствии с распределением Пуассона. Зависимая переменная должна быть неотрицательной.

Функции связи между таргетом и объясняющими переменными предполагается логарифмической: $$g(\eta) = \ln(\eta) \Rightarrow \hat {y} = e^{x{^T}\beta + {\beta_{0}}}.$$

Модель оценивается методом максимального правдоподобия, функция логарифма правдоподобия с учетом штрафа регуляризации эластичной сети имеет вид: $$\max_{\beta,\beta_0} \frac{1}{N} \sum_{i=1}^{N} \Big( y_i(x_{i}^{T}\beta + \beta_0) - e^{x{^T_i}\beta + {\beta_0}} \Big)- \lambda \Big( \alpha||\beta||_1 + \dfrac {1} {2}(1 - \alpha)||\beta||^2_2 \Big),$$

где

  • $\lambda$ – параметр, отвечающий за силу регуляризации. $\lambda\in\mathbb{R}^{+}$;
  • $\alpha$ – параметр, отвечающий за распределение штрафов регуляризации между нормой 1 ($\ell_1$) и нормой 2 ($\ell_2$). $\alpha\in[0,1]$;
  • $||\beta||{_1}$ – штраф регуляризации $\ell_1$ (LASSO). $||\beta||{_1} = \sum{^p_{k=1}} |\beta{_k}|$;
  • $||\beta||{_2}$ – штраф регуляризации $\ell_2$ (Ridge). $||\beta||{_2} = \sum{^p_{k=1}} \beta{^2_k}$.

Тогда соответствующая метрика Deviance имеет вид: $$D = -2 \sum_{i=1}^{N} \big( y_i \text{ln}(y_i / \hat {y}_i) - (y_i - \hat {y}_i) \big).$$

Вывод функции правдоподобия для GLM с распределением Пуассона (без регуляризации)

NB: Ниже $\lambda$ не имеет отношения к вышеупомянутому одноименному параметру регуляризации.

Напомним, что функция вероятности для распределения Пуассона имеет вид: $$p(k;\lambda) = \frac{\lambda^k e^{-\lambda}}{k!}, \hspace{10pt} \lambda\in\mathbb{R}^{+}.$$ Также, для распределения Пуассона справедливо, что: $$\mathbb{E}\left[k\right] = Var(k) = \lambda.$$ Тогда, для оценивания коэффициентов нашей модели необходимо максимизировать правдоподобие (совместную условную вероятность при имеющихся данных), что данные имеют распределение Пуассона: $$p(y_1,\dots,y_n|x_1,\dots,x_n;\beta_0,\beta) = \prod_{i=1}^{N}\frac{e^{y_i(x_i{^T}\beta + {\beta_{0}})} e^{-e^{x_i{^T}\beta + {\beta_{0}}}}}{y_i!} = L(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n).$$ Для упрощения задачи оптимизации перейдем к логарифму правдоподобия: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(y_i(x_i{^T}\beta + {\beta_{0}}) -e^{x_i{^T}\beta + {\beta_{0}}}-\ln(y_i!)\right).$$ Поскольку величина $\ln(y_i!)$ не зависит от выбора параметров, можно упростить задачу: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(y_i(x_i{^T}\beta + {\beta_{0}}) -e^{x_i{^T}\beta + {\beta_{0}}}\right).$$ Далее решается задача оптимизации для определения параметров модели (пакеты реализуют оптимизацию численными методами): $$\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta_0} = 0,\\\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta} = 0.$$ Обычно минимизируется отрицательное правдоподобие, которое является выпуклой функцией.

Вывод метрики Deviance для GLM с распределением Пуассона

Метрика Deviance представляет собой отношение правдоподобия между двумя моделями: рассматриваемой моделью и "идеальной" моделью, в которая бы идеально предсказывала бы зависимую переменную. $$Deviance = 2(\ell_{ideal} - \ell_{model})$$

В качестве такой "идеальной модели" может использоваться сама зависимая переменная. Тогда, логарифм правдоподобия "идеальной модели" для GLM с распределением Пуассона имеет вид: $$\ell_{ideal} = \sum_{i=1}^{N}\left(y_i \ln(y_i) -y_i-\ln(y_i!)\right).$$

Из приведенного выше вывода правдоподобия для рассматриваемой модели, мы можем записать, обозначив $\hat{y}_i = e^{x{^T}\beta + {\beta_{0}}}$: $$\ell_{model} = \sum_{i=1}^{N}\left(y_i \ln(\hat{y}_i) -\hat{y}_i-\ln(y_i!)\right).$$

Тогда получаем, $$Deviance = 2\sum_{i=1}^{N}\left(y_i \ln(y_i) -y_i - y_i \ln(\hat{y}_i) +\hat{y}_i\right) = -2\sum_{i=1}^{N}\left(y_i \ln(y_i/\hat{y}_i) - (y_i -\hat{y}_i)\right).$$

GLM с гамма-распределением

GLM с гамма-распределением используется для моделирования положительной непрерывной зависимой переменной, когда ее условная дисперсия увеличивается вместе со средним значением, но коэффициент вариации зависимой переменной предполагается постоянным.

Обычно GLM с гамма-распределением используются с логарифмической или обратной функциями связи: $$g(\eta) = \ln(\eta);\hspace{20pt}g(\eta) = \frac{1}{\eta}.$$

Модель оценивается методом максимального правдоподобия, функция логарифма правдоподобия (для обратной функции связи) с учетом штрафа регуляризации эластичной сети имеет вид: $$\max_{\beta,\beta_0} - \frac{1}{N} \sum_{i=1}^{N} \frac{y_i}{x{^T_i}\beta + \beta_0} + \text{ln} \big( x{^T_i}\beta + \beta_0 \big ) - \lambda \Big( \alpha||\beta||_1 + \dfrac {1} {2}(1 - \alpha)||\beta||^2_2 \Big),$$

где

  • $\lambda$ – параметр, отвечающий за силу регуляризации. $\lambda\in\mathbb{R}^{+}$;
  • $\alpha$ – параметр, отвечающий за распределение штрафов регуляризации между нормой 1 ($\ell_1$) и нормой 2 ($\ell_2$). $\alpha\in[0,1]$;
  • $||\beta||{_1}$ – штраф регуляризации $\ell_1$ (LASSO). $||\beta||{_1} = \sum{^p_{k=1}} |\beta{_k}|$;
  • $||\beta||{_2}$ – штраф регуляризации $\ell_2$ (Ridge). $||\beta||{_2} = \sum{^p_{k=1}} \beta{^2_k}$.

Соответствующая метрика Deviance имеет вид: $$D = 2 \sum_{i=1}^{N} - \text{ln} \bigg (\dfrac {y_i} {\hat {y}_i} \bigg) + \dfrac {(y_i - \hat{y}_i)} {\hat {y}_i}.$$

Вывод функции правдоподобия для GLM с Гамма-распределением и логарифмической функцией связи (без регуляризации)

NB: Ниже $\alpha$ и $\lambda$ не имеют отношения к вышеупомянутым одноименным параметрам регуляризации.

Напомним, что функция плотности вероятности для Гамма-распределения имеет вид: $$f(x;\alpha,\lambda) = \begin{cases} \frac{\alpha^\lambda}{\Gamma(\lambda)}x^{\lambda-1}e^{-\alpha x},&x \ge 0\\ 0,& x <0 \end{cases},\hspace{20pt} \Gamma(x) = \int_0^{\infty}x^{\lambda-1}e^{-x}dx. $$ Также, для Гамма-распределения справедливо, что: $$ \mathbb{E}[x] = \frac{\lambda}{\alpha},\hspace{15pt} Var(x) = \frac{\lambda}{\alpha^2}.$$

Для оценивания GLM удобно параметризовать данное распределение иначе: $$ \mathbb{E}[x] = \frac{\lambda}{\alpha} = \mu,\hspace{15pt} Var(x) = \frac{\lambda}{\alpha^2} = \frac{\mu^2}{\lambda}.$$ Тогда, $$f(x;\mu,\lambda) = \frac{1}{\Gamma(\lambda)}\alpha^\lambda x^{\lambda-1}e^{-\alpha x} = \frac{1}{\Gamma(\lambda)}\left(\frac{\lambda}{\mu}\right)^\lambda x^{\lambda-1}e^{-\left(\frac{\lambda}{\mu}\right)x} = \frac{1}{x\cdot\Gamma(\lambda)}\left(\frac{\lambda x}{\mu}\right)^\lambda e^{-\frac{\lambda x}{\mu}},\hspace{15pt} x\ge 0.$$

Также напомним, что для GLM с функцией связи $\ln(x)$ предсказание (оценка математического ожидания) имеет вид $\hat {y} = e^{x{^T}\beta + {\beta_{0}}}$.

Тогда, для оценивания коэффициентов нашей модели необходимо максимизировать правдоподобие (совместную условную вероятность при имеющихся данных), что данные имеют Гамма-распределение, в предположении известного параметра $\lambda$: $$p(y_1,\dots,y_n|x_1,\dots,x_n;\beta_0,\beta) = \prod_{i=1}^{N}\frac{1}{y_i\cdot\Gamma(\lambda)}\left(\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)^\lambda e^{-\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}} = L(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n).$$ Для упрощения задачи оптимизации перейдем к логарифму правдоподобия: \begin{align} \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) &= \sum_{i=1}^{N}\left(-\ln(\Gamma(\lambda))-\ln(y_i)+ \lambda\ln\left(\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)\\ &= \sum_{i=1}^{N}\left(-\ln(\Gamma(\lambda))-\ln(y_i) + \lambda\ln\left(\lambda\right)+\lambda\ln\left(y_i\right)-\lambda\left({x_i{^T}\beta + {\beta_{0}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)\\ &= \sum_{i=1}^{N}\left((\lambda-1)\ln\left(y_i\right)-\ln(\Gamma(\lambda)) + \lambda\ln\left(\lambda\right)-\lambda\left({x_i{^T}\beta + {\beta_{0}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right). \end{align} Тогда обозначим слагаемые, не зависящие от параметров модели ($\beta_0$, $\beta$), как некоторую константу $C$, получаем: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(-\lambda\left(\left({x_i{^T}\beta + {\beta_{0}}}\right) +\frac{y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right) + C(y_i, \lambda)\right),$$ Таким образом, требуется максимизировать $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = -\sum_{i=1}^{N}\left(\left({x_i{^T}\beta + {\beta_{0}}}\right) + \frac{y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right),$$ Далее решается задача оптимизации для определения параметров модели (пакеты реализуют оптимизацию численными методами): $$\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta_0} = 0,\\\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta} = 0.$$ Обычно минимизируется отрицательное правдоподобие, которое является выпуклой функцией.

Вывод метрики Deviance для GLM с Гамма-распределением

Метрика Deviance представляет собой отношение правдоподобия между двумя моделями: рассматриваемой моделью и "идеальной" моделью, в которая бы идеально предсказывала бы зависимую переменную. $$Deviance = 2(\ell_{ideal} - \ell_{model})$$

В качестве такой "идеальной модели" может использоваться сама зависимая переменная. Тогда, логарифм правдоподобия "идеальной модели" для GLM с распределением Пуассона имеет вид: $$\ell_{ideal} = -\sum_{i=1}^{N}\left(1 + \ln(y_i)\right).$$

Из приведенного выше вывода правдоподобия для рассматриваемой модели, мы можем записать, обозначив $\hat{y}_i = e^{x{^T}\beta + {\beta_{0}}}$: $$\ell_{model} = -\sum_{i=1}^{N}\left(\frac{y_i}{\hat{y}_i} + \ln(\hat{y}_i)\right).$$

Тогда получаем, $$Deviance = 2\sum_{i=1}^{N}\left(- 1 - \ln(y_i) + \frac{y_i}{\hat{y}_i} + \ln(\hat{y}_i)\right) = 2\sum_{i=1}^{N}\left(-\ln\left(\frac{y_i}{\hat{y}_i}\right)+\frac{y_i-\hat{y}_i}{\hat{y}_i} \right).$$

Дополнительная литература по GLM

Установка H2O на Google Colaboratory и инициализация

In [31]:
!apt-get install default-jre
Reading package lists... Done
Building dependency tree       
Reading state information... Done
default-jre is already the newest version (2:1.11-68ubuntu1~18.04.1).
default-jre set to manually installed.
The following package was automatically installed and is no longer required:
  libnvidia-common-440
Use 'apt autoremove' to remove it.
0 upgraded, 0 newly installed, 0 to remove and 35 not upgraded.
In [12]:
!java -version
openjdk version "1.8.0_265"
OpenJDK Runtime Environment (build 1.8.0_265-b01)
OpenJDK 64-Bit Server VM (build 25.265-b01, mixed mode)
In [30]:
pip install h2o
Collecting h2o
  Downloading h2o-3.30.1.1.tar.gz (129.3 MB)
     |████████████████████████████████| 129.3 MB 36 kB/s  eta 0:00:01
Requirement already satisfied: requests in /home/egor/anaconda3/envs/py_3.8.5_j/lib/python3.8/site-packages (from h2o) (2.24.0)
Collecting tabulate
  Downloading tabulate-0.8.7-py3-none-any.whl (24 kB)
Requirement already satisfied: future in /home/egor/anaconda3/envs/py_3.8.5_j/lib/python3.8/site-packages (from h2o) (0.18.2)
Collecting colorama>=0.3.8
  Downloading colorama-0.4.3-py2.py3-none-any.whl (15 kB)
Requirement already satisfied: idna<3,>=2.5 in /home/egor/anaconda3/envs/py_3.8.5_j/lib/python3.8/site-packages (from requests->h2o) (2.10)
Requirement already satisfied: certifi>=2017.4.17 in /home/egor/anaconda3/envs/py_3.8.5_j/lib/python3.8/site-packages (from requests->h2o) (2020.6.20)
Requirement already satisfied: chardet<4,>=3.0.2 in /home/egor/anaconda3/envs/py_3.8.5_j/lib/python3.8/site-packages (from requests->h2o) (3.0.4)
Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /home/egor/anaconda3/envs/py_3.8.5_j/lib/python3.8/site-packages (from requests->h2o) (1.25.10)
Building wheels for collected packages: h2o
  Building wheel for h2o (setup.py) ... done
  Created wheel for h2o: filename=h2o-3.30.1.1-py2.py3-none-any.whl size=129358603 sha256=57226983b329bc677b3cc2c06877667f0a50da8345cbb17c163313e3eeee4c2f
  Stored in directory: /home/egor/.cache/pip/wheels/bd/35/7e/d7b4a455bdbc7e30afce970b6ce471282015c4aa13415ea670
Successfully built h2o
Installing collected packages: tabulate, colorama, h2o
Successfully installed colorama-0.4.3 h2o-3.30.1.1 tabulate-0.8.7
Note: you may need to restart the kernel to use updated packages.
In [31]:
import h2o
from h2o.estimators.glm import H2OGeneralizedLinearEstimator
h2o.init()
Checking whether there is an H2O instance running at http://localhost:54321 ..... not found.
Attempting to start a local H2O server...
  Java Version: openjdk version "10.0.2" 2018-07-17; OpenJDK Runtime Environment Zulu10.3+5 (build 10.0.2+13); OpenJDK 64-Bit Server VM Zulu10.3+5 (build 10.0.2+13, mixed mode)
  Starting server from /home/egor/anaconda3/envs/py_3.8.5_j/lib/python3.8/site-packages/h2o/backend/bin/h2o.jar
  Ice root: /tmp/tmp45mnkhk_
  JVM stdout: /tmp/tmp45mnkhk_/h2o_egor_started_from_python.out
  JVM stderr: /tmp/tmp45mnkhk_/h2o_egor_started_from_python.err
  Server is running at http://127.0.0.1:54321
Connecting to H2O server at http://127.0.0.1:54321 ... successful.
H2O_cluster_uptime: 01 secs
H2O_cluster_timezone: Europe/Moscow
H2O_data_parsing_timezone: UTC
H2O_cluster_version: 3.30.1.1
H2O_cluster_version_age: 5 days
H2O_cluster_name: H2O_from_python_egor_6aieos
H2O_cluster_total_nodes: 1
H2O_cluster_free_memory: 5.857 Gb
H2O_cluster_total_cores: 6
H2O_cluster_allowed_cores: 6
H2O_cluster_status: accepting new members, healthy
H2O_connection_url: http://127.0.0.1:54321
H2O_connection_proxy: {"http": null, "https": null}
H2O_internal_security: False
H2O_API_Extensions: Amazon S3, XGBoost, Algos, AutoML, Core V3, TargetEncoder, Core V4
Python_version: 3.8.5 final

Построение GLM для частоты страховых случаев

In [32]:
# Преобразование в H2O-Frame

h2o_train_c = h2o.H2OFrame(pd.concat([x_train_c, y_train_c], axis=1))
h2o_valid_c = h2o.H2OFrame(pd.concat([x_valid_c, y_valid_c], axis=1))
h2o_test_c = h2o.H2OFrame(pd.concat([x_test_c, y_test_c], axis=1))
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
In [33]:
# Инициализируем и обучим GLM модель c кросс-валидацией

glm_poisson = H2OGeneralizedLinearEstimator(family = "poisson", link = "Log", nfolds=5)
glm_poisson.train(y="ClaimsCount", x = h2o_train_c.names[1:-1], training_frame = h2o_train_c, validation_frame = h2o_valid_c)
glm Model Build progress: |███████████████████████████████████████████████| 100%
In [34]:
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных

glm_poisson.summary()
GLM Model: summary
family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame
0 poisson log Elastic Net (alpha = 0.5, lambda = 9.667E-5 ) 21 19 3 Key_Frame__upload_8f8b855e97db5670f288e3fd8ba141d9.hex
Out[34]:

In [35]:
# Метрики качества модели - по всем данным и на кросс-валидации

glm_poisson.cross_validation_metrics_summary().as_data_frame()
Out[35]:
mean sd cv_1_valid cv_2_valid cv_3_valid cv_4_valid cv_5_valid
0 mae 0.39099473 0.0027418558 0.38808832 0.39292654 0.3933419 0.38791674 0.3927002
1 mean_residual_deviance 1.0612507 0.010642719 1.0458891 1.0696331 1.0717574 1.0556002 1.0633734
2 mse 0.5886037 0.011173075 0.57106256 0.5964719 0.596917 0.5838192 0.5947478
3 null_deviance 17458.027 250.67941 17174.996 17692.336 17652.361 17200.85 17569.592
4 r2 0.008045993 0.0028944598 0.004718711 0.010126316 0.007916696 0.0058058705 0.01166237
5 residual_deviance 17109.195 199.60555 16902.613 17293.828 17328.174 16923.38 17097.98
6 rmse 0.7671773 0.0073083583 0.7556868 0.7723159 0.77260405 0.76408064 0.7711989
7 rmsle 0.35901806 0.0018053672 0.35685304 0.3603561 0.36082768 0.35734305 0.35971043
In [36]:
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)

glm_poisson._model_json['output']['coefficients_table'].as_data_frame()
Out[36]:
names coefficients standardized_coefficients
0 Intercept -2.565659 -1.575867
1 Gender -0.163708 -0.079349
2 MariStat -0.114494 -0.041245
3 HasKmLimit -0.439625 -0.137290
4 BonusMalus 0.012823 0.196838
5 OutUseNb 0.081984 0.057037
6 RiskArea 0.010119 0.022421
7 driver_minage_m 0.000209 0.003928
8 driver_minage_f 0.012369 0.198390
9 driver_minage_m_2 0.000000 0.000000
10 driver_minage_f_2 -0.000091 -0.113901
11 VehUsage_Private -0.190581 -0.090123
12 VehUsage_Private+trip to office 0.000000 0.000000
13 VehUsage_Professional 0.320065 0.105780
14 VehUsage_Professional run 0.491840 0.065987
15 SocioCateg_CSP1 -0.079219 -0.012003
16 SocioCateg_CSP2 -0.205329 -0.034207
17 SocioCateg_CSP3 0.242698 0.025177
18 SocioCateg_CSP4 0.050770 0.012696
19 SocioCateg_CSP5 -0.003721 -0.001769
20 SocioCateg_CSP6 0.060435 0.024822
21 SocioCateg_CSP7 -0.495631 -0.005237
In [37]:
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации

pmodels = {}
pmodels['overall'] = glm_poisson.coef_norm()
for x in range(len(glm_poisson.cross_validation_models())):
    pmodels[x] = glm_poisson.cross_validation_models()[x].coef_norm()
coef = pd.DataFrame.from_dict(pmodels).round(5)
coef['overall'] = abs(coef['overall'])
coef.sort_values('overall',ascending=False)
Out[37]:
overall 0 1 2 3 4
Intercept 1.57587 -1.56929 -1.58006 -1.58025 -1.57437 -1.57666
driver_minage_f 0.19839 0.19686 0.26079 0.09732 0.09520 0.12461
BonusMalus 0.19684 0.20471 0.18426 0.19930 0.20297 0.19012
HasKmLimit 0.13729 -0.12588 -0.13843 -0.15366 -0.13748 -0.13194
driver_minage_f_2 0.11390 -0.10674 -0.17633 -0.02598 -0.01771 -0.05999
VehUsage_Professional 0.10578 0.10780 0.11021 0.10484 0.10246 0.10419
VehUsage_Private 0.09012 -0.10194 -0.08889 -0.07295 -0.10099 -0.08568
Gender 0.07935 -0.07373 -0.10003 -0.06666 -0.05018 -0.05652
VehUsage_Professional run 0.06599 0.07869 0.05949 0.06946 0.06185 0.05947
OutUseNb 0.05704 0.04973 0.05546 0.05601 0.06296 0.06115
MariStat 0.04125 -0.03105 -0.04759 -0.03853 -0.04047 -0.05028
SocioCateg_CSP2 0.03421 -0.05051 -0.02785 -0.03595 -0.03389 -0.02617
SocioCateg_CSP3 0.02518 0.02744 0.02781 0.02495 0.02337 0.02200
SocioCateg_CSP6 0.02482 0.03324 0.01906 0.00703 0.02664 0.03006
RiskArea 0.02242 0.02870 0.02384 0.01873 0.01966 0.02106
SocioCateg_CSP4 0.01270 0.01618 0.00924 0.00918 0.01658 0.01069
SocioCateg_CSP1 0.01200 -0.00999 -0.01034 -0.01074 -0.01373 -0.01552
SocioCateg_CSP7 0.00524 -0.00425 -0.00040 -0.00435 -0.03253 -0.00364
driver_minage_m 0.00393 0.01993 0.00000 0.00000 0.01317 0.00014
SocioCateg_CSP5 0.00177 0.00000 0.00000 0.00000 -0.00240 -0.00824
driver_minage_m_2 0.00000 0.00000 -0.00943 0.00175 0.00000 0.00000
VehUsage_Private+trip to office 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
In [38]:
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок

c_train_pred = glm_poisson.predict(h2o_train_c).as_data_frame()
c_valid_pred = glm_poisson.predict(h2o_valid_c).as_data_frame()
c_test_pred = glm_poisson.predict(h2o_test_c).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
In [39]:
# Сохранение обученной модели

model_glm_poisson = h2o.save_model(model=glm_poisson, path="data/", force=True)
In [40]:
model_glm_poisson
Out[40]:
'/home/egor/Projects/GeekBrains/EDA/7_lesson/data/GLM_model_python_1597496210036_1'

Построение GLM для среднего убытка

In [41]:
# Преобразование в H2O-Frame

h2o_train_ac = h2o.H2OFrame(pd.concat([x_train_ac, y_train_ac], axis=1))
h2o_valid_ac = h2o.H2OFrame(pd.concat([x_valid_ac, y_valid_ac], axis=1))
h2o_test_ac = h2o.H2OFrame(pd.concat([x_test_ac, y_test_ac], axis=1))
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
In [87]:
h2o_train_ac
driver_minexp Gender MariStat HasKmLimit BonusMalus OutUseNb RiskArea driver_minage_m driver_minage_f driver_minage_m_2 driver_minage_f_2 VehUsage_Private VehUsage_Private+trip to office VehUsage_Professional VehUsage_Professional run SocioCateg_CSP1 SocioCateg_CSP2 SocioCateg_CSP3 SocioCateg_CSP4 SocioCateg_CSP5 SocioCateg_CSP6 SocioCateg_CSP7 AvgClaim
52 0 0 0 50 0 7 38 18 1444 324 0 1 0 0 0 0 0 0 1 0 0 42.7292
52 0 0 0 50 2 7 56 18 3136 324 0 1 0 0 0 0 0 1 0 0 0 520.785
52 0 1 0 50 0 7 59 18 3481 324 1 0 0 0 0 0 0 0 0 1 0 927.481
52 1 0 1 50 0 6 18 65 324 4225 1 0 0 0 0 0 0 0 0 1 0 372.609
52 1 0 0 50 0 6 18 49 324 2401 0 0 1 0 0 0 0 0 1 0 01568.17
52 1 0 0 95 1 11 18 63 324 3969 0 1 0 0 0 0 0 0 1 0 0 9.77356
52 1 1 1 76 1 10 18 33 324 1089 1 0 0 0 0 0 0 0 1 0 01973.62
45 0 1 0 90 0 6 25 18 625 324 0 1 0 0 0 0 0 0 1 0 03086.02
52 0 0 0 50 0 6 44 18 1936 324 0 1 0 0 0 0 0 0 1 0 0 709.305
52 1 1 0 85 0 10 18 43 324 1849 0 1 0 0 0 0 0 0 1 0 0 390.589
Out[87]:

In [42]:
# Инициализируем и обучим GLM модель c кросс-валидацией

glm_gamma = H2OGeneralizedLinearEstimator(family = "gamma", link = "Log", nfolds=5)
glm_gamma.train(y="AvgClaim", x = h2o_train_ac.names[1:-1], training_frame = h2o_train_ac, validation_frame = h2o_valid_ac)
glm Model Build progress: |███████████████████████████████████████████████| 100%
In [43]:
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных

glm_gamma.summary()
GLM Model: summary
family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame
0 gamma log Elastic Net (alpha = 0.5, lambda = 2.602E-4 ) 21 19 3 Key_Frame__upload_9d98cd5d45dd69fb647507022407470b.hex
Out[43]:

In [44]:
# Метрики качества модели - по всем данным и на кросс-валидации

glm_gamma.cross_validation_metrics_summary().as_data_frame()
Out[44]:
mean sd cv_1_valid cv_2_valid cv_3_valid cv_4_valid cv_5_valid
0 mae 1192.6599 78.264984 1074.1938 1152.9423 1257.785 1226.5664 1251.8121
1 mean_residual_deviance 2.1104307 0.18365581 1.8613946 1.9738325 2.2194502 2.2034786 2.293997
2 mse 1.5812803E7 1.9704548E7 4330197.5 6515766.5 1.1213453E7 6235834.5 5.0768764E7
3 null_deviance 3177.9753 260.33975 2851.8447 2983.9773 3475.8562 3206.393 3371.8057
4 r2 -0.004497921 0.0065900427 -0.008070656 -0.001890941 0.0026067737 -0.01408022 -0.0010545618
5 residual_deviance 3209.7488 265.3541 2870.2705 2992.33 3460.1228 3294.2004 3431.8198
6 rmse 3520.9094 2066.3977 2080.9126 2552.5999 3348.6494 2497.1653 7125.22
7 rmsle 1.7666154 0.046983425 1.7337893 1.7003165 1.8003181 1.7929715 1.8056816
In [45]:
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)

glm_gamma._model_json['output']['coefficients_table'].as_data_frame()
Out[45]:
names coefficients standardized_coefficients
0 Intercept 7.310145 7.017288
1 Gender -0.314623 -0.152506
2 MariStat 0.299338 0.110533
3 HasKmLimit -0.066859 -0.016827
4 BonusMalus -0.002345 -0.039579
5 OutUseNb -0.021085 -0.015749
6 RiskArea 0.030396 0.067861
7 driver_minage_m -0.021733 -0.397791
8 driver_minage_f 0.012161 0.187589
9 driver_minage_m_2 0.000179 0.268589
10 driver_minage_f_2 -0.000154 -0.183893
11 VehUsage_Private 0.000000 0.000000
12 VehUsage_Private+trip to office 0.024089 0.012019
13 VehUsage_Professional -0.021973 -0.008023
14 VehUsage_Professional run -0.290266 -0.044476
15 SocioCateg_CSP1 -0.129253 -0.020325
16 SocioCateg_CSP2 -0.218260 -0.037997
17 SocioCateg_CSP3 -0.069526 -0.006869
18 SocioCateg_CSP4 -0.019099 -0.005290
19 SocioCateg_CSP5 0.000000 0.000000
20 SocioCateg_CSP6 0.092551 0.035955
21 SocioCateg_CSP7 -1.358740 -0.015578
In [46]:
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации

pmodels = {}
pmodels['overall'] = glm_gamma.coef_norm()
for x in range(len(glm_gamma.cross_validation_models())):
    pmodels[x] = glm_gamma.cross_validation_models()[x].coef_norm()
coef_ac = pd.DataFrame.from_dict(pmodels).round(5)
coef_ac['overall'] = abs(coef['overall'])
coef_ac.sort_values('overall',ascending=False)
Out[46]:
overall 0 1 2 3 4
Intercept 1.57587 7.03779 7.01743 7.00762 7.00037 7.01081
driver_minage_f 0.19839 0.33056 0.04862 -0.00838 0.22504 0.30055
BonusMalus 0.19684 -0.05567 -0.04954 -0.05031 -0.05661 0.01878
HasKmLimit 0.13729 -0.03149 -0.00742 -0.00620 -0.00957 -0.02741
driver_minage_f_2 0.11390 -0.29734 -0.07856 -0.01410 -0.21851 -0.28017
VehUsage_Professional 0.10578 0.00000 -0.00766 -0.01850 -0.03836 0.00868
VehUsage_Private 0.09012 -0.01263 0.00000 0.00000 0.00000 -0.00015
Gender 0.07935 -0.19715 -0.08128 -0.01080 -0.20057 -0.26181
VehUsage_Professional run 0.06599 -0.05041 -0.03016 -0.06644 -0.04736 -0.03271
OutUseNb 0.05704 -0.02159 -0.00919 -0.00540 -0.03373 -0.00313
MariStat 0.04125 0.13655 0.13818 0.11761 0.12983 0.01112
SocioCateg_CSP2 0.03421 -0.05301 -0.06124 -0.04068 -0.02647 -0.02822
SocioCateg_CSP3 0.02518 -0.00496 -0.01733 -0.01289 0.00000 -0.01090
SocioCateg_CSP6 0.02482 0.04000 0.00000 0.02141 0.01318 0.04691
RiskArea 0.02242 0.07030 0.05127 0.07033 0.06968 0.08519
SocioCateg_CSP4 0.01270 -0.00153 -0.05073 -0.01323 -0.00637 0.00405
SocioCateg_CSP1 0.01200 -0.02541 -0.02537 -0.03700 -0.00783 -0.02360
SocioCateg_CSP7 0.00524 -0.01695 0.00000 -0.01552 -0.01741 -0.01882
driver_minage_m 0.00393 -0.47229 -0.17808 -0.01291 -0.40483 -0.95331
SocioCateg_CSP5 0.00177 0.00000 -0.05960 0.00000 0.00000 0.00000
driver_minage_m_2 0.00000 0.34617 0.09684 -0.01871 0.21061 0.73959
VehUsage_Private+trip to office 0.00000 0.00825 0.02511 0.00041 0.00294 0.00000
In [47]:
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок

ac_train_pred = glm_gamma.predict(h2o_train_ac).as_data_frame()
ac_valid_pred = glm_gamma.predict(h2o_valid_ac).as_data_frame()
ac_test_pred = glm_gamma.predict(h2o_test_ac).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
In [48]:
# Сохранение обученной модели

model_glm_gamma = h2o.save_model(model=glm_gamma, path="data/", force=True)
In [49]:
model_glm_gamma
Out[49]:
'/home/egor/Projects/GeekBrains/EDA/7_lesson/data/GLM_model_python_1597496210036_2'

Использование GLM моделей

In [50]:
h2o_df = h2o.H2OFrame(df_freq[col_features])
Parse progress: |█████████████████████████████████████████████████████████| 100%
In [51]:
df['CountPredicted'] = glm_poisson.predict(h2o_df).as_data_frame()
df['AvgClaimPredicted'] = glm_gamma.predict(h2o_df).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
In [53]:
df['BurningCost'] = df.CountPredicted * df.AvgClaimPredicted
df[['CountPredicted', 'AvgClaimPredicted', 'BurningCost']].head()
Out[53]:
CountPredicted AvgClaimPredicted BurningCost
index
145780 NaN NaN NaN
145781 NaN NaN NaN
145782 NaN NaN NaN
145783 NaN NaN NaN
145784 NaN NaN NaN

* Домашнее задание: GLM для прогнозирования наступления страхового случая

In [54]:
# Разбиение датасета на train/val/test

df_ind = data.get_pd(col_features+['ClaimInd'])

x_train_ind, x_test_ind, y_train_ind, y_test_ind = train_test_split(df_ind[col_features], df_ind.ClaimInd, test_size=0.3, random_state=1)
x_valid_ind, x_test_ind, y_valid_ind, y_test_ind = train_test_split(x_test_ind, y_test_ind, test_size=0.5, random_state=1)
In [55]:
# Преобразование в H2O-Frame

h2o_train_c = h2o.H2OFrame(pd.concat([x_train_ind, y_train_ind], axis=1))
h2o_valid_c = h2o.H2OFrame(pd.concat([x_valid_ind, y_valid_ind], axis=1))
h2o_test_c = h2o.H2OFrame(pd.concat([x_test_ind, y_test_ind], axis=1))
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
In [76]:
# Преобразуем целевую переменную ClaimInd в категориальную при помощи метода asfactor во всех наборах данных

h2o_train_c['ClaimInd'].asfactor()
h2o_valid_c['ClaimInd'].asfactor()
h2o_test_c['ClaimInd'].asfactor()
Out[76]:
['driver_minexp',
 'Gender',
 'MariStat',
 'HasKmLimit',
 'BonusMalus',
 'OutUseNb',
 'RiskArea',
 'driver_minage_m',
 'driver_minage_f',
 'driver_minage_m_2',
 'driver_minage_f_2',
 'VehUsage_Private',
 'VehUsage_Private+trip to office',
 'VehUsage_Professional',
 'VehUsage_Professional run',
 'SocioCateg_CSP1',
 'SocioCateg_CSP2',
 'SocioCateg_CSP3',
 'SocioCateg_CSP4',
 'SocioCateg_CSP5',
 'SocioCateg_CSP6',
 'SocioCateg_CSP7',
 'ClaimInd']
In [89]:
h2o_train_c['ClaimInd']
ClaimInd
1
1
0
0
0
0
0
0
0
0
Out[89]:

In [79]:
# Инициализируем и обучим GLM модель c кросс-валидацией
In [78]:
glm_binomial = H2OGeneralizedLinearEstimator(family = "binomial", link = "Logit", nfolds=5)
In [91]:
# Инициализируем и обучим GLM модель c кросс-валидацией

glm_binomial.train(y="ClaimInd", x = h2o_train_c.names[1:-1], training_frame = h2o_train_c, validation_frame = h2o_valid_c)
glm Model Build progress: |███████████████████████████████████████████████| 100%
In [93]:
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных

glm_binomial.summary()
GLM Model: summary
family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame
0 binomial logit Elastic Net (alpha = 0.5, lambda = 2.368E-5 ) 21 21 3 Key_Frame__upload_a3f4d416bedf8fd89e04fe1259d4f41.hex
Out[93]:

In [94]:
# Метрики качества модели - по всем данным и на кросс-валидации

glm_binomial.cross_validation_metrics_summary().as_data_frame()
Out[94]:
mean sd cv_1_valid cv_2_valid cv_3_valid cv_4_valid cv_5_valid
0 accuracy 0.5004415 0.08289997 0.44190806 0.47953686 0.40942773 0.5654677 0.605867
1 auc 0.56346345 0.009767907 0.5569546 0.57068455 0.5497133 0.5725396 0.56742525
2 aucpr 0.113697775 0.0037565897 0.11407095 0.119015664 0.10844074 0.11309471 0.113866806
3 err 0.4995585 0.08289997 0.55809194 0.5204631 0.59057224 0.43453228 0.39413297
4 err_count 8052.2 1324.3737 8997.0 8406.0 9484.0 7019.0 6355.0
5 f0point5 0.13302433 0.0064273607 0.12966333 0.13869794 0.12400163 0.13339314 0.13936567
6 f1 0.18759094 0.0065006316 0.18571816 0.19667432 0.17901662 0.1862029 0.19034272
7 f2 0.31904125 0.015040828 0.3271475 0.33793104 0.32177755 0.30822968 0.30012053
8 lift_top_group 1.4246209 0.14822257 1.224569 1.4448918 1.4580715 1.6312454 1.3643265
9 logloss 0.31126234 0.0062752496 0.31403515 0.3188719 0.30980614 0.30180112 0.31179735
10 max_per_class_error 0.5409678 0.06405939 0.58166975 0.53886294 0.6192799 0.4526244 0.5124021
11 mcc 0.0571638 0.011409857 0.049557984 0.06553044 0.040806193 0.06632031 0.06360409
12 mean_per_class_accuracy 0.5480664 0.009996953 0.541419 0.5547651 0.5338816 0.5573253 0.55294096
13 mean_per_class_error 0.4519336 0.009996953 0.458581 0.44523486 0.46611837 0.44267473 0.44705904
14 mse 0.085404895 0.0022214192 0.08631573 0.088193156 0.08474387 0.08214188 0.08562982
15 null_deviance 10103.041 203.34572 10179.154 10380.477 9990.729 9839.38 10125.466
16 pr_auc 0.113697775 0.0037565897 0.11407095 0.119015664 0.10844074 0.11309471 0.113866806
17 precision 0.111442305 0.0061621494 0.107943185 0.11591754 0.10291629 0.11218217 0.11825234
18 r2 0.003908301 0.0011621623 0.003313827 0.0046493188 0.0022391728 0.0051963087 0.0041428762
19 recall 0.60698354 0.08549763 0.66450775 0.6483932 0.6870432 0.5473756 0.4875979
20 residual_deviance 10036.101 204.47351 10125.122 10300.199 9950.354 9749.987 10054.841
21 rmse 0.2922213 0.0038084106 0.29379538 0.29697332 0.291108 0.28660405 0.29262573
22 specificity 0.48914927 0.10041951 0.41833025 0.46113706 0.38072008 0.5672749 0.618284
In [95]:
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)

glm_binomial._model_json['output']['coefficients_table'].as_data_frame()
Out[95]:
names coefficients standardized_coefficients
0 Intercept -2.475348e+00 -2.279640
1 Gender -6.416985e-02 -0.031103
2 MariStat -6.469029e-02 -0.023304
3 HasKmLimit -3.698560e-01 -0.115502
4 BonusMalus 6.792484e-03 0.104267
5 OutUseNb 6.145252e-02 0.042753
6 RiskArea 9.240660e-03 0.020474
7 driver_minage_m -3.727295e-03 -0.070187
8 driver_minage_f -1.168549e-03 -0.018743
9 driver_minage_m_2 -6.145260e-06 -0.009595
10 driver_minage_f_2 -9.652992e-08 -0.000121
11 VehUsage_Private -1.437545e-01 -0.067979
12 VehUsage_Private+trip to office -1.656392e-02 -0.008276
13 VehUsage_Professional 2.260091e-01 0.074695
14 VehUsage_Professional run 2.965853e-01 0.039791
15 SocioCateg_CSP1 -1.073627e-01 -0.016267
16 SocioCateg_CSP2 -1.462515e-01 -0.024365
17 SocioCateg_CSP3 1.043504e-01 0.010825
18 SocioCateg_CSP4 2.983183e-02 0.007460
19 SocioCateg_CSP5 -2.803307e-02 -0.013326
20 SocioCateg_CSP6 4.853385e-02 0.019934
21 SocioCateg_CSP7 -2.029539e-01 -0.002144
In [96]:
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации

pmodels = {}
pmodels['overall'] = glm_binomial.coef_norm()
for x in range(len(glm_binomial.cross_validation_models())):
    pmodels[x] = glm_binomial.cross_validation_models()[x].coef_norm()
coef_ac = pd.DataFrame.from_dict(pmodels).round(5)
coef_ac['overall'] = abs(coef['overall'])
coef_ac.sort_values('overall',ascending=False)
Out[96]:
overall 0 1 2 3 4
Intercept 1.57587 -2.28394 -2.28940 -2.27901 -2.26686 -2.28048
driver_minage_f 0.19839 -0.01369 -0.01345 -0.08917 -0.01683 0.00000
BonusMalus 0.19684 0.11242 0.10418 0.10726 0.09011 0.10216
HasKmLimit 0.13729 -0.12437 -0.10507 -0.12688 -0.10837 -0.11304
driver_minage_f_2 0.11390 0.02622 -0.00099 0.04655 -0.01538 -0.02174
VehUsage_Professional 0.10578 0.07059 0.07935 0.07858 0.07359 0.08750
VehUsage_Private 0.09012 -0.04901 -0.06051 -0.07380 -0.07220 -0.06452
Gender 0.07935 -0.07189 -0.03813 -0.01354 -0.01978 -0.05341
VehUsage_Professional run 0.06599 0.04448 0.04789 0.03218 0.03362 0.04745
OutUseNb 0.05704 0.04975 0.04220 0.03977 0.03916 0.04315
MariStat 0.04125 -0.02557 -0.02110 -0.03040 -0.02404 -0.01776
SocioCateg_CSP2 0.03421 -0.02134 -0.03077 -0.02779 -0.01912 -0.01725
SocioCateg_CSP3 0.02518 0.00576 0.00822 0.02285 0.01346 0.00927
SocioCateg_CSP6 0.02482 0.00000 0.01309 0.03802 0.02743 0.03027
RiskArea 0.02242 0.01686 0.01640 0.03170 0.02688 0.01079
SocioCateg_CSP4 0.01270 0.01207 -0.00191 0.02384 0.00809 0.00629
SocioCateg_CSP1 0.01200 -0.01809 -0.02563 -0.01074 -0.01472 -0.00614
SocioCateg_CSP7 0.00524 -0.00260 -0.03202 0.00343 0.00187 -0.00267
driver_minage_m 0.00393 -0.17758 -0.06282 -0.13569 -0.04993 -0.16287
SocioCateg_CSP5 0.00177 -0.00046 -0.01264 -0.00924 -0.01366 -0.00957
driver_minage_m_2 0.00000 0.08494 -0.02226 0.05205 -0.03821 0.07889
VehUsage_Private+trip to office 0.00000 -0.00535 -0.00718 0.00000 -0.00671 0.00000
In [105]:
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок

с_train_pred = glm_binomial.predict(h2o_train_c).as_data_frame()
с_valid_pred = glm_binomial.predict(h2o_valid_c).as_data_frame()
с_test_pred = glm_binomial.predict(h2o_test_c).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%

Сохраним полученную модель

In [103]:
model_glm_binomial = h2o.save_model(model=glm_binomial, path="data/", force=True)
In [112]:
с_train_pred.iloc[:, 0]
Out[112]:
0        0
1        0
2        1
3        1
4        0
        ..
80603    0
80604    0
80605    0
80606    0
80607    1
Name: predict, Length: 80608, dtype: int64
In [119]:
from sklearn.metrics import accuracy_score, f1_score, confusion_matrix
In [120]:
# Выведем импортированные выше метрики классификации для обучающей, валидационной и тестовой выборок

cm_train_binominal = confusion_matrix(y_train_ind ,с_train_pred.iloc[:, 0])
cm_valid_binominal = confusion_matrix(y_valid_ind ,с_valid_pred.iloc[:, 0])
cm_test_binominal = confusion_matrix(y_test_ind ,с_test_pred.iloc[:, 0])

cm_valid_binominal, cm_train_binominal, cm_test_binominal
Out[120]:
(array([[9240, 6424],
        [ 784,  825]]),
 array([[42814, 30159],
        [ 3741,  3894]]),
 array([[9168, 6481],
        [ 801,  824]]))

Как видно из матриц выше, у нас большое количество FN значений по всем данным.

In [123]:
f1_train_binominal = f1_score(y_train_ind ,с_train_pred.iloc[:, 0], average='weighted')
f1_valid_binominal = f1_score(y_valid_ind ,с_valid_pred.iloc[:, 0], average='weighted')
f1_test_binominal = f1_score(y_test_ind ,с_test_pred.iloc[:, 0], average='weighted')

print(f'train - {f1_train_binominal}, valid - {f1_valid_binominal}, test - {f1_test_binominal}')
train - 0.6662250021448088, valid - 0.6697403995228527, test - 0.6657756605965353
In [124]:
accuracy_score

ac_train_binominal = accuracy_score(y_train_ind ,с_train_pred.iloc[:, 0])
ac_valid_binominal = accuracy_score(y_valid_ind ,с_valid_pred.iloc[:, 0])
ac_test_binominal = accuracy_score(y_test_ind ,с_test_pred.iloc[:, 0])

print(f'train - {ac_train_binominal}, valid - {ac_valid_binominal}, test - {ac_test_binominal}')
train - 0.579446208813021, valid - 0.5827013257685405, test - 0.5784415885145305

Метрики показывают, что модель получилось явно слабой и близка просто к угадыванию 50% на 50%. Думаю, что ситуацию можно исправить сгенерировать новые фичи. Так же выборка является не сбалансированной, что тоже плохо сказалось на модели.